#### Functions are invoked by the main replication file
extract <- function(y){### To compute balance statistics
  c(mean.tr=y$mean.Tr,
  mean.ctrl=y$mean.Co,
  sd.diff=y$sdiff.pooled/10,
#  eQQMed=y$qqsummary$mediandiff,
#    eQQMean=y$qqsummary$meandiff,
#    eQQMax=y$qqsummary$maxdiff,
  var.ratio=y$var.ratio,
  t.test=y$p.value,
  ks=if(is.null(y$ks$ks$statistic)){NA}else{y$ks$ks$statistic},
  ks.p.boot=if(is.null(y$ks$ks.boot.pvalue)){NA}else{y$ks$ks.boot.pvalue})}


stack.tlm <- function(x){#x is a summary(lm()) object 
	## Stacking function for experimental regressions
	## Uses only the intercept plus variables with TRUE in the name drops contro
	if(class(x)!="summary.lm"){cat("Object of wrong class! \n")}else{
	 		select <- grep("Intercept|TRUE|PP|TT|cue",row.names(x$coef)) 	
 core<-stack(data.frame(round(rbind(x$coef[select,"Estimate"],
 		x$coef[select,"Std. Error"],x$coef[select,"Pr(>|t|)"]),3)))
 core$ind<-as.character(core$ind)
 names(core)[1]<-"lm"
 core$type <- rep(1:3,nrow(core)/3)
 core$ind <- gsub("\\.","",core$ind,perl=T)
 core$ind <- gsub("TT_somecue","TREAT",core$ind,perl=T)
 core$ind <- gsub("^I.*\\d\\d\\.TRUE","PARTY",core$ind,perl=T)
 core$ind <- gsub("P13TRUE","PT",core$ind,perl=T)
 core$ind <- gsub("P45TRUE","PSDB",core$ind,perl=T)
 core$ind <- gsub("PP13","PT",core$ind,perl=T)
 core$ind <- gsub("PP45","PSDB",core$ind,perl=T)
 core$ind <- gsub("^XIntercept","INTERCEPT",core$ind,perl=T)
 core<-rbind(core,c(as.numeric(nrow(x$coef)>length(select)),"zControls",1))
 core<-rbind(core,c(length(x$residuals),"zN",2))
 core<-rbind(core,c(round(x$r.squared,3),"zR2",3))
 core[,1]<-as.numeric(core[,1])  #to keep numbers as numbers  
 return(core)}}
 
test.diff.effects <- function(x,pty){ 
	### Used for head on comparison between pt cue and psdb cues
	### X is a bootstrapped object (matrix with 1000 estimates)
	#function uses bootstraped estimates to test difference
	pticuenteraction <- paste("PP",pty,":T2ptcue",sep="")
	psdbcueinteraction <- paste("PP",pty,":T3psdbcue",sep="")
	the.diff<- x[,"T2ptcue"]+x[,pticuenteraction]-#pt cue effect on ptistas
			  (x[,"T3psdbcue"]+x[,psdbcueinteraction])#psdb cue effect on ptistas
	pe.diff <-mean(the.diff)
	se.diff <-sd(the.diff)
	p.diff <- sum(sign(the.diff)==sign(pe.diff))/length(the.diff)
	cat("Reporting PT cue effect - PSDB cue effect\n")
	out <- c(diff=pe.diff,se=se.diff,pval=p.diff)
		}	
		
exp.sometreat <- function(qq,partyvar="PP",bt=T,the.data=d){ 
				#qq is the name of the experimental question
				#Can be applied to both pilot and lapop experimental data
				#Runs a "single" treatment condition (some treatment)
				#Reports regressions with and without controls for PT and PSDB
				#Linear probability models without controls, and logit with controls
				#Uses PP as the definition of party
	vars <- paste(c("Idade","Sexo","iA","iK")[is.element(
				  c("Idade","Sexo","iA","iK"),names(the.data))])
  	lcl.data <- subset(the.data,select=c("TT","PP",qq,vars))
	cat("\n--------------Estimation Details ------------\n")
		cat("Missing observations by variable:\n")
	print(as.matrix(apply(is.na(lcl.data),2,sum)))
	lcl.data$agree.pt <- lcl.data[,qq]==13
	lcl.data$agree.psdb <- lcl.data[,qq]==45
	if(partyvar!="PP"){lcl.data$PP <- the.data[,partyvar]}	
	cat("\tRunning Agreement ~ Party + Party * SomeCue\n")
	cat("\tParty levels:",levels(lcl.data$PP),"\n");
	cat("\tQuestion:",qq,"\n")
	controls <-   paste(vars,collapse="+")
  	form.pt <- as.formula(paste("agree.pt~PP*TT",controls,sep="+"))
	form.psdb <- as.formula(paste("agree.psdb~PP*TT",controls,sep="+"))
	cat("\tControls:",controls,"\n") #Select control variables that exist in data set
	cat("\nN Party Identification by Treatment Condition")	
	print(table(lcl.data$PP))
	cat("\n\tEstimating...\n")
	reg.pt <- lm(agree.pt~ PP*TT,data=lcl.data) #PP is just PT, PSDB and others
  	reg.psdb <- lm(agree.psdb~ PP*TT,data=lcl.data)
  	reg.ptC <-glm(formula=form.pt,data=lcl.data,family = binomial("logit"),x=T)
	#reg.ptC <- zelig(form.pt,model="logit",data=tmp,cite=F) 		  
	reg.psdbC <- glm(formula=form.psdb,data=lcl.data,family = binomial("logit"),x=T)
  	### Bootstrap SE's
  	if(bt==T){
  	my.stat.13 <- function(ddd,i){ #function for boostrapping PT 
  		ddd <- ddd[i,]
  		reg <- lm(agree.pt~ PP*TT,data=ddd)
  		coefficients(reg)
  	}
  	my.stat.45 <- function(ddd,i){ #function for boostrapping PSDB
  		ddd <- ddd[i,]
  		reg <- lm(agree.psdb~ PP*TT,data=ddd)
  		coefficients(reg)
  	}
  	require(boot)
  	cat("\n\tBootstrapping linear probability model (takes time)...\n")
  	my.boot13 <- boot(lcl.data, my.stat.13, 1000)
  	se.boot13 <- apply(my.boot13[[2]],2,sd)
  	se.boot13 <- my.boot13[[2]]
  	colnames(se.boot13) <- names(coef(reg.pt))
  	my.boot45 <- boot(lcl.data, my.stat.45, 1000)
  	se.boot45 <- apply(my.boot45[[2]],2,sd)
  	se.boot45 <- my.boot45[[2]]
  	colnames(se.boot45) <- names(coef(reg.psdb))
  }else{se.boot13<-se.boot45<-NULL}#
  if(bt==T){#Needed because Zelig no longer runs well...
  		### Bootstrap SE's
  		my.stat.13C <- function(ddd,i){ #function for boostrapping PT (no controls)
  		ddd <- ddd[i,]
  		reg <- glm(formula=form.pt,data=ddd,family = binomial("logit"),x=T)
  		coefficients(reg)
  		}
  		my.stat.45C <- function(ddd,i){ #function for boostrapping PSDB (no controls)
  		ddd <- ddd[i,]
  		reg <- glm(formula=form.psdb,data=ddd,family = binomial("logit"),x=T)
  		coefficients(reg)
  		}
  		cat("\tBootstrapping logit models (takes time)...\n")
  		require(boot)
  		my.boot13C <- boot(lcl.data, my.stat.13C, 1000)
  		#se.boot13C <- apply(my.boot13C[[2]],2,sd)
  		se.boot13C <- my.boot13C[[2]]
  		colnames(se.boot13C) <- names(coef(reg.ptC))
   		my.boot45C <- boot(lcl.data, my.stat.45C, 1000)
  		#se.boot45C <- apply(my.boot45C[[2]],2,sd)
  		se.boot45C <- my.boot45C[[2]]
  		colnames(se.boot45C) <- names(coef(reg.psdbC))
  	}else{se.boot13C<-se.boot45C<-NULL}#
  #Create table of linear probability model
	tmp.table <- merge(stack.tlm(summary(reg.pt)),
				  stack.tlm(summary(reg.psdb)),
				  by=c('ind','type'),suffixes=c("PT","PSDB"))
 	tmp.table$ind <- ifelse(tmp.table$type==2,"SE",
 					ifelse(tmp.table$type==3,"p-value",tmp.table$ind))
 	tmp.table$ind[c(nrow(tmp.table)-1,nrow(tmp.table))] <- c("zN","zR$^2$")
 	tables <- tmp.table[,-2]
  to.return<-list(tables=tables,reg13ctrl=reg.ptC,reg13=reg.pt,
  				reg45ctrl=reg.psdbC,reg45=reg.psdb,
  				reg13boot=se.boot13,reg45boot=se.boot45,
  				reg13ctrlboot=se.boot13C,reg45ctrlboot=se.boot45C)
  cat("---------------------Done--------------------\n\n")
 return(to.return)
}


exp.eachtreat <- function(qq,partyvar="PP",bt=T,the.data=d){ 
	require(car)
					#qq is the name of the experimental question
					#Applied to  FB experimental data
					#Runs all treatment condition (eacb treatment)
					#Currently reports regressions with and without controls 
					#for agreement with PT and PSDB
					#Linear probability models without controls, logits with controls
					#Uses PP (definition of party, just PT, PSDB, and everybodyelse)	  
	vars <- paste(c("Idade","Sexo","iA")[is.element(
				  c("Idade","Sexo","iA"),names(the.data))]) 
				   
 	lcl.data <- subset(the.data,select=c("T","PP",qq,vars))
 	#tmp.dk <-tmp
 	if(partyvar!="PP"){lcl.data$PP <- the.data[,partyvar]}
	#lcl.data <- lcl.data[which(d[,qq]!=77),] #eliminate the don't knows
	lcl.data$agree.pt <- lcl.data[,qq]==13
	lcl.data$agree.psdb <- lcl.data[,qq]==45
	cat("\n--------------Estimation Details ------------\n")
	cat("Estmating: Agreement ~ Party + Party * EachCue\n")
	cat("\tParty levels:",levels(lcl.data$PP),"\n");
	cat("\tQuestion:",qq,"\n")
	controls <-   paste(vars,collapse="+")
  	form.pt <- as.formula(paste("agree.pt~PP*T",controls,sep="+"))
	form.psdb <- as.formula(paste("agree.psdb~PP*T",controls,sep="+"))
	cat("\tControls:",controls,"\n") 
	cat("\n\tN Party Identification by Treatment Condition")
	print(table(lcl.data$PP,lcl.data$T))
  	cat("\n\tEstimating...\n")
	reg.pt <- lm(agree.pt~ PP*T,data=lcl.data) 	
	#reg.ptCz <- zelig(formula=form.pt,model="logit",data=lcl.data,cite=F) 
	reg.ptC <-glm(formula=form.pt,data=lcl.data,family = binomial("logit"),x=T)#x=T need for erer ma
	reg.psdb <- lm(agree.psdb~ PP*T,data=lcl.data)	
	#reg.psdbC <-  zelig(form.psdb,model="logit",data=lcl.data,cite=F)
	reg.psdbC <-  glm(form.psdb,data=lcl.data,family = binomial("logit"),x=T)
  	if(bt==T){
  		### Bootstrap SE's
  		my.stat.13 <- function(ddd,i){ #function for boostrapping PT (no controls)
  		ddd <- ddd[i,]
  		reg <- lm(agree.pt~ PP*T,data=ddd)
  		coefficients(reg)
  		}
  		my.stat.45 <- function(ddd,i){ #function for boostrapping PSDB (no controls)
  		ddd <- ddd[i,]
  		reg <- lm(agree.psdb~ PP*T,data=ddd)
  		coefficients(reg)
  		}
  		cat("\tBootstrapping linear probability model (takes time)...\n")
  		require(boot)
  		my.boot13 <- boot(lcl.data, my.stat.13, 1000)
  		#se.boot13 <- apply(my.boot13[[2]],2,sd)#this is just SEs
  		se.boot13 <- my.boot13[[2]]
  		colnames(se.boot13) <- names(coef(reg.pt))
   		my.boot45 <- boot(lcl.data, my.stat.45, 1000)
  		#se.boot45 <- apply(my.boot45[[2]],2,sd)#this is just SEs
  		se.boot45 <- my.boot45[[2]]
  		colnames(se.boot45) <- names(coef(reg.psdb))
  	}else{se.boot13<-se.boot45<-NULL}#
  	if(bt==T){#Needed because Zelig no longer runs well...
  		### Bootstrap SE's
  		my.stat.13C <- function(ddd,i){ #function for boostrapping PT (no controls)
  		ddd <- ddd[i,]
  		reg <- glm(formula=form.pt,data=ddd,family = binomial("logit"),x=T)
  		coefficients(reg)
  		}
  		my.stat.45C <- function(ddd,i){ #function for boostrapping PSDB (no controls)
  		ddd <- ddd[i,]
  		reg <- glm(formula=form.psdb,data=ddd,family = binomial("logit"),x=T)
  		coefficients(reg)
  		}
  		cat("\tBootstrapping logit models (takes time)...\n")
  		require(boot)
  		my.boot13C <- boot(lcl.data, my.stat.13C, 1000)
  		#se.boot13C <- apply(my.boot13C[[2]],2,sd)
  		se.boot13C <- my.boot13C[[2]]
  		colnames(se.boot13C) <- names(coef(reg.ptC))
   		my.boot45C <- boot(lcl.data, my.stat.45C, 1000)
  		#se.boot45C <- apply(my.boot45C[[2]],2,sd)
  		se.boot45C <- my.boot45C[[2]]
  		colnames(se.boot45C) <- names(coef(reg.psdbC))
  	}else{se.boot13C<-se.boot45C<-NULL}#
	#Assemble regression table
	tab <- merge(stack.tlm(summary(reg.pt)),
				  stack.tlm(summary(reg.psdb)),
				  by=c('ind','type'),suffixes=c("PT","PSDB"))
	tab$ind <- ifelse(tab$type==2,"SE",ifelse(tab$type==3,"p-value",tab$ind))
	tab$ind[c(nrow(tab)-1,nrow(tab))] <- c("zN","zR$^2$")
	tables <- tab[,-2] 
 	to.return<-list(reg13ctrl=reg.ptC,reg13=reg.pt,
 					reg45ctrl=reg.psdbC,reg45=reg.psdb,
 					reg13boot=se.boot13,reg45boot=se.boot45,
 					reg13ctrlboot=se.boot13C,reg45ctrlboot=se.boot45C)
 cat("---------------------Done--------------------\n\n")
 return(to.return)
}


effects.cue <- function(xx,boot.se=NULL){
	#Computes treatment effects and SE from lm, glmm, and boostraped objects
	#x is a list of regressions
	#boot.se is an object with bootstrapped SE results
	interactionspe <-c("PP45:T2ptcue","PP45:T3psdbcue","PP45:T4dblcue",
				   "PP13:T2ptcue","PP13:T3psdbcue","PP13:T4dblcue","PP66:T4dblcue")
	basepe <- c("T2ptcue","T3psdbcue","T4dblcue","T2ptcue",
				"T3psdbcue","T4dblcue","T4dblcue")
	if(class(xx[[1]])[1]=="lm"){ #linear models (wihtout controls)
	cat("Treatment effects based on linear probability models\n")
	#Point estimates for agreement with the party
	pes <- lapply(xx,function(x){
		t1<-coef(x)[interactionspe];	
		t0<-coef(x)[basepe];
		tt<-t1+t0;
		names(tt)<-c(interactionspe);
		return(tt)}) 
	ses <- lapply(xx,function(x){
		var1<-diag(vcov(x))[interactionspe];
		var0<-diag(vcov(x))[basepe]	
		cov01<-c(vcov(x)[interactionspe[1],basepe[1]],
			 vcov(x)[interactionspe[2],basepe[2]],
			 vcov(x)[interactionspe[3],basepe[3]],
			 vcov(x)[interactionspe[4],basepe[4]],
			 vcov(x)[interactionspe[5],basepe[5]],
			 vcov(x)[interactionspe[6],basepe[6]],
			 vcov(x)[interactionspe[7],basepe[7]])
		ss<- sqrt(var1+var0+2*cov01)
		return(ss)})
	pes0 <- lapply(xx,function(x){ #no interaction term (non partisans)
		out <- coef(x)[c("T2ptcue","T3psdbcue","T4dblcue")]
		names(out) <- c('PP0:T2ptcue','PP0:T3psdbcue','PP0:T4dblcue')
		return(out)
	  })
	ses0 <-  lapply(xx,function(x){ #no interaction term (non partisans)
		out <- sqrt(diag(vcov(x))[c("T2ptcue","T3psdbcue","T4dblcue")])
		names(out) <- c('PP0:T2ptcue','PP0:T3psdbcue','PP0:T4dblcue')
		return(out)
	  })
	baseline <- t(sapply(xx,function(x){#get baselines for bar graphs
			c("PP0"=sum(coef(x)[1]),
			  "PP45"=sum(coef(x)[1],coef(x)['PP45']),
			  "PP13"=sum(coef(x)[1],coef(x)['PP13']),
			  "PP66"=sum(coef(x)[1],coef(x)['PP66']))}))
	pes <- cbind(do.call('rbind',pes0),do.call('rbind',pes))
	ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))
	if(is.null(boot.se)==F){ ##bootstrapped linear models
		cat("Using bootstrapped SE\n")  
		ses <- lapply(boot.se,function(x){apply(x[,interactionspe]+x[,basepe],2,sd)})
		ses0 <-  lapply(boot.se,function(x){
						out<-apply(x[,c("T2ptcue","T3psdbcue","T4dblcue")],2,sd)
						names(out)<-c('PP0:T2ptcue','PP0:T3psdbcue','PP0:T4dblcue');
						return(out)})
		ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))
	}	
}## end linear

	if(class(xx[[1]])[1]=="glm"){ ### logit models (with controls)	
	cat("Treatment effects based on logit models:
		 Point Estimates computed using 'erer' package, which corrects for,
		 corrects for binary independent variables, and averages marginal effects
		 over all values of other variables instead of holding them at average")
	require(erer)#loads maBina
	#Point estimates for agreement with the party
	pes <- lapply(xx,function(x){
		me <- maBina(x, x.mean = T, rev.dum = T, digits = 5)$out  
		t1<-me[interactionspe,"effect"];	
		t0<-me[basepe,"effect"];
		tt<-t1+t0;
		names(tt)<-c(interactionspe);
		return(tt)}) 
	ses <- lapply(pes,function(x)is.na(x))#placeholdder
	pes0 <- lapply(xx,function(x){ #no interaction term (non partisans)
		me <- maBina(x, x.mean = T, rev.dum = T, digits = 5)$out  
		out <- me[c("T2ptcue","T3psdbcue","T4dblcue"),"effect"]
		names(out) <- c('PP0:T2ptcue','PP0:T3psdbcue','PP0:T4dblcue')
		return(out)
	  })
	ses0 <-  lapply(xx,function(x){ #no interaction term (non partisans)
		me <- maBina(x, x.mean = T, rev.dum = T, digits = 5)$out  
		out<-me[c("T2ptcue","T3psdbcue","T4dblcue"),"error"]
		names(out) <- c('PP0:T2ptcue','PP0:T3psdbcue','PP0:T4dblcue')
		return(out)
	  })
	pes <- cbind(do.call('rbind',pes0),do.call('rbind',pes))
	ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))		
	## Baseline predicted values 
	## No treatment, non-identifier, women, average age and iA
	baseline <- t(sapply(xx,function(x){#get baseline values for plot
			m.Idade <- mean(x$model$Idade)
			m.iA <- mean(x$model$iA)
			m.Sexo <- factor("Masculino")
		x_base <- data.frame(
				PP=factor(levels(x$model$PP)),
				T="_ctrl",
				Idade=m.Idade,Sexo=m.Sexo,iA=m.iA)
		out<- predict(x,newdata=x_base,type = "response",se.fit=T)$fit
		names(out) <- paste("PP",levels(x$model$PP),sep="")
	    return(out)}))
	baseline <- baseline[,c(1,3,2,4)]##for backwards compatibility
	if(is.null(boot.se)==F){
		cat("Using bootstrapped SE\n")  
		partype <- c("PP45","PP45","PP45","PP13","PP13","PP13","PP66")
		ses <- lapply(boot.se,function(x){apply(inv.logit(
												x[,interactionspe]+
												x[,basepe]+
												x[,partype]+
												x[,"(Intercept)"]),2,sd)})
		ses0 <-  lapply(boot.se,function(x){
						out<-apply(inv.logit(
									x[,c("T2ptcue","T3psdbcue","T4dblcue")]+
									x[,"(Intercept)"]),2,sd)			
						names(out)<-c('PP0:T2ptcue','PP0:T3psdbcue','PP0:T4dblcue');
						return(out)})
		ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))
	}else{
	cat("\nFor SE, declare boot.se=allTQs.PARTY.ctrl.bootse\n")}
	}#end GLS	
	#collect and finish
	rownames(pes) <- rownames(ses) <- rownames(baseline) <- Qs
	effects <- list(pes=pes,ses=ses,baseline=baseline) 
	return(effects)
}


effects.pooled <- function(xx,boot.se=NULL){
	### Takes in either a list of lm or glm objects 
	### Nead to deal with bootstrapping still...
	### Compute effects of pooled treatment
	interactionspe <-c("PP45:TT_somecue","PP13:TT_somecue","PP66:TT_somecue")
	basepe <- c("TT_somecue")	
	# Regular Point estimates for agreement with the PT
	if(class(xx[[1]])[1]=="lm"){ #linear models (wihtout controls)
	cat("Pooled-treatment effects based on linear probability models\n")
	pes <- lapply(xx,function(x){
		t1<-coef(x)[interactionspe];	
		t0<-coef(x)[basepe];
		tt<-t1+t0;
		names(tt)<-interactionspe;return(tt)}) 
	ses <- lapply(xx,function(x){
		var1<-diag(vcov(x))[interactionspe];
		var0<-diag(vcov(x))[basepe]	
		cov01<-	c(vcov(x)[interactionspe[1],basepe],
			 vcov(x)[interactionspe[2],basepe],
			 vcov(x)[interactionspe[3],basepe])
		ss<- sqrt(var1+var0+2*cov01)
		return(ss)})
	pes0 <- lapply(xx,function(x){
  			t1<-coef(x)[basepe];	
			names(t1)<-"PP0:TT_somecue";return(t1)})
	ses0 <- lapply(xx,function(x){
			se<-sqrt(diag(vcov(x))[basepe])	
			names(se)<-"PP0:TT_somecue"
			return(se)})
	baseline <- t(sapply(xx,function(x){#get baselines for bar graphs
			c("PP0"=sum(coef(x)[1]),
			  "PP45"=sum(coef(x)[1],coef(x)['PP45']),
			  "PP13"=sum(coef(x)[1],coef(x)['PP13']),
			  "PP66"=sum(coef(x)[1],coef(x)['PP66'])
			  )}))
	pes <- cbind(do.call('rbind',pes0),do.call('rbind',pes))
	ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))
	if(is.null(boot.se)==F){ ##bootstrapped models
		cat("Using bootstrapped SE\n")  
		ses <- lapply(boot.se,function(x){apply(x[,interactionspe]+x[,basepe],2,sd)})
		ses0 <-  lapply(boot.se,function(x){out<-sd(x[,basepe])
									names(out)<-"PP0:TT_somecue";
										return(out)})
		ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))
			}#end lm.boostrap
	}#end of lm	
	if(class(xx[[1]])[1]=="glm"){ ### logit models (with controls)	
	cat("Treatment effects based on logit models:
		 Point Estimates computed using 'erer' package, which corrects for,
		 corrects for binary independent variables, and averages marginal effects
		 over all values of other variables instead of holding them at average")
	require(erer)#loads maBina
	#Point estimates for agreement with the party
	pes <- lapply(xx,function(x){
		me <- maBina(x, x.mean = T, rev.dum = T, digits = 5)$out  
		t1<-me[interactionspe,"effect"];	
		t0<-me[basepe,"effect"];
		tt<-t1+t0;
		names(tt)<-c(interactionspe);
		return(tt)}) 
	ses <- lapply(pes,function(x)is.na(x))#placeholdder
	pes0 <- lapply(xx,function(x){ #no interaction term (non partisans)
		me <- maBina(x, x.mean = T, rev.dum = T, digits = 5)$out  
		out <- me[c("TT_somecue"),"effect"]
		names(out) <- c('PP0:TT_somecue')
		return(out)
	  })
	ses0 <-  lapply(xx,function(x){ #no interaction term (non partisans)
		me <- maBina(x, x.mean = T, rev.dum = T, digits = 5)$out  
		out<-me[c("TT_somecue"),"error"]
		names(out) <- c('PP0:TT_somecue')
		return(out)
	  })
	pes <- cbind(do.call('rbind',pes0),do.call('rbind',pes))
	ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))		
	## Baseline predicted values 
	## No treatment, non-identifier, women, average age and iA
	baseline <- t(sapply(xx,function(x){#get baseline values for plot
			m.Idade <- mean(x$model$Idade)
			m.iA <- mean(x$model$iA)
			m.Sexo <- factor("Masculino")
		x_base <- data.frame(
				PP=factor(levels(x$model$PP)),
				TT="_ctrl",
				Idade=m.Idade,Sexo=m.Sexo,iA=m.iA)
		out<- predict(x,newdata=x_base,type = "response",se.fit=T)$fit
		names(out) <- paste("PP",levels(x$model$PP),sep="")
	    return(out)}))
	baseline <- baseline[,c(1,3,2,4)]##for backwards compatibility
	if(is.null(boot.se)==F){#bootstrapped SE for glm
		cat("Using bootstrapped SE\n")  
		partype <- c("PP45","PP13","PP66")
		ses <- lapply(boot.se,function(x){apply(inv.logit(
												x[,interactionspe]+
												x[,basepe]+
												x[,partype]+
												x[,"(Intercept)"]),2,sd)})
		ses0 <-  lapply(boot.se,function(x){
						out<-sd(inv.logit(
									x[,c("TT_somecue")]+
									x[,"(Intercept)"]))			
						names(out)<-c('PP0:TT_somecue');
						return(out)})
		ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))
	}else{
	cat("\n\nFor SE, declare boot.se=allTQs.PARTY.ctrl.bootse\n")}
	}#end GLS	
#	pes <- cbind(do.call('rbind',pes0),do.call('rbind',pes))
#	ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))
#	baseline <- do.call('rbind',baseline)
	rownames(pes) <- rownames(ses) <-  rownames(baseline) <-names(xx)
	effects <- list(pes=pes,ses=ses,baseline=baseline)
	return(effects)
}


effects.pooled.other <- function(xx,boot.se=NULL){
	### Takes in either a list of lm or glm objects 
	### Compute effects of pooled treatment
	### For other parties only
	interactionspe <-c("PP15:TT_somecue","PP43:TT_somecue","PP66:TT_somecue")
	basepe <- c("TT_somecue")	
	# Regular Point estimates for agreement with the PT
	if(class(xx[[1]])[1]=="lm"){ #linear models (wihtout controls)
	cat("Pooled-treatment effects based on linear probability models\n")
	pes <- lapply(xx,function(x){
		t1<-coef(x)[interactionspe];	
		t0<-coef(x)[basepe];
		tt<-t1+t0;
		names(tt)<-interactionspe;return(tt)}) 
	ses <- lapply(xx,function(x){
		var1<-diag(vcov(x))[interactionspe];
		var0<-diag(vcov(x))[basepe]	
		cov01<-	c(vcov(x)[interactionspe[1],basepe],
			 vcov(x)[interactionspe[2],basepe],
			 vcov(x)[interactionspe[3],basepe])
		ss<- sqrt(var1+var0+2*cov01)
		return(ss)})
	pes0 <- lapply(xx,function(x){
  			t1<-coef(x)[basepe];	
			names(t1)<-"PP0:TT_somecue";return(t1)})
	ses0 <- lapply(xx,function(x){
			se<-sqrt(diag(vcov(x))[basepe])	
			names(se)<-"PP0:TT_somecue"
			return(se)})
	baseline <- t(sapply(xx,function(x){#get baselines for bar graphs
			c("PP0"=sum(coef(x)[1]),
			  "PP15"=sum(coef(x)[1],coef(x)['PP15']),
			  "PP43"=sum(coef(x)[1],coef(x)['PP43']),
			  "PP66"=sum(coef(x)[1],coef(x)['PP66'])
			  )}))
	if(is.null(boot.se)==F){ ##bootstrapped models
		cat("Using bootstrapped SE\n")  
		ses <- lapply(boot.se,function(x){apply(x[,interactionspe]+x[,basepe],2,sd)})
		ses0 <-  lapply(boot.se,function(x){out<-sd(x[,basepe])
									names(out)<-"PP0:TT_somecue";
										return(out)})
			}
	pes <- cbind(do.call('rbind',pes0),do.call('rbind',pes))
	ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))
	}#end of lm	
	if(class(xx[[1]])[1]=="glm"){ ### logit models (with controls)	
	cat("Pooled-treatment effects based on logit models:
		 Point Estimates computed using 'erer' package, which corrects for,
		 corrects for binary independent variables, and averages marginal effects
		 over all values of other variables instead of holding them at average")
	require(erer)#loads maBina
	#Point estimates for agreement with the party
	pes <- lapply(xx,function(x){
		me <- maBina(x, x.mean = T, rev.dum = T, digits = 5)$out  
		t1<-me[interactionspe,"effect"];	
		t0<-me[basepe,"effect"];
		tt<-t1+t0;
		names(tt)<-c(interactionspe);
		return(tt)}) 
	ses <- lapply(pes,function(x)is.na(x))#placeholdder
	pes0 <- lapply(xx,function(x){ #no interaction term (non partisans)
		me <- maBina(x, x.mean = T, rev.dum = T, digits = 5)$out  
		out <- me[c("TT_somecue"),"effect"]
		names(out) <- c('PP0:TT_somecue')
		return(out)
	  })
	ses0 <-  lapply(xx,function(x){ #no interaction term (non partisans)
		me <- maBina(x, x.mean = T, rev.dum = T, digits = 5)$out  
		out<-me[c("TT_somecue"),"error"]
		names(out) <- c('PP0:TT_somecue')
		return(out)
	  })
	pes <- cbind(do.call('rbind',pes0),do.call('rbind',pes))
	ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))		
	## Baseline predicted values 
	## No treatment, non-identifier, women, average age and iA
	baseline <- t(sapply(xx,function(x){#get baseline values for plot
			m.Idade <- mean(x$model$Idade)
			m.iA <- mean(x$model$iA)
			m.Sexo <- factor("Masculino")
		x_base <- data.frame(
				PP=factor(levels(x$model$PP)),
				TT="_ctrl",
				Idade=m.Idade,Sexo=m.Sexo,iA=m.iA)
		out<- predict(x,newdata=x_base,type = "response",se.fit=T)$fit
		names(out) <- paste("PP",levels(x$model$PP),sep="")
	    return(out)}))
	baseline <- baseline[,c(1,3,4,6)]##for backwards compatibility
	if(is.null(boot.se)==F){#bootstrapped SE for glm
		cat("Using bootstrapped SE\n")  
		partype <- c("PP15","PP43","PP66")
		ses <- lapply(boot.se,function(x){apply(inv.logit(
												x[,interactionspe]+
												x[,basepe]+
												x[,partype]+
												x[,"(Intercept)"]),2,sd)})
		ses0 <-  lapply(boot.se,function(x){
						out<-sd(inv.logit(
									x[,c("TT_somecue")]+
									x[,"(Intercept)"]))			
						names(out)<-c('PP0:TT_somecue');
						return(out)})
		ses <- cbind(do.call('rbind',ses0),do.call('rbind',ses))
	}else{
	cat("\n\nFor SE, declare boot.se=allTQs.PARTY.ctrl.bootse\n")}
	}#end of glm
	rownames(pes) <- rownames(ses) <-  rownames(baseline) <-Qs
	effects <- list(pes=pes,ses=ses,baseline=baseline)
	return(effects)
}


plot.effects.pooled <- function(xx,leg=T,leg2=T,ommit.others=F,port=F){
	#Takes the ouput of effects.pooled
	pes <- xx[["pes"]] 
	ses <- xx[["ses"]]
	offsets <- c(.07,.21,-.21,-.07)
	pchs <- c(23,21,21,21) 
	cols <- c(gray(0.5),1,1,gray(0.5))
	ltys <- c(1,1,1)
	bgs <- c(gray(0.5),1,"white",gray(0.5)) 
	if(ommit.others==T){
		others <- grep("66",colnames(pes))
		pes <- pes[,-others]
		ses <- ses[,-others]
		offsets <- c(0,.14,-.14)
		pchs <- pchs[-others]
		cols <- cols[-others]
		ltys <-ltys[-others]
		bgs <- bgs[-others]
	}
	par(mar=c(7,8,.7,1))
	ylims <- c(0.5,0.5+nrow(pes))
	xlims <- c(-0.5,0.5)
	plot(xlims,ylims,type="n",ylab="",xlab="",yaxt="n")
	my.polygon<-function(x,col){#background
	polygon(x=c(xlims,xlims[2:1]),
			y=c(rep(x-0.5,2),rep(x+0.5,2)),col=gray(col),border=NA)
	}
	mapply(my.polygon,1:nrow(pes),
			rep(c(0.8, 0.7),nrow(pes))[1:nrow(pes)])
	for(i in 1:nrow(pes)){#lines
		cilow <- pes[i,]-qnorm(.975)*ses[i,]
		cihigh <- pes[i,]+qnorm(.975)*ses[i,]
		segments(x0=cilow,x1=cihigh,y0=i+offsets,y1=i+offsets,
				lwd=2,col=cols,lty=ltys)
		points(pes[i,],i+offsets,pch=pchs,bg=bgs,col=cols,cex=2)}		
	abline(v=0,lty=2)
	if(port==F){
	if(leg){axis(side=2,at=1:5,tick=F,
			labels=c("Venezuela\nin Mercosul","Offshore Oil\n(Pre-salt)",
					 "Gvt. Financing\nPvt. Companies","Investment\nProtection",
					 "Minimum\nWage"),las=2,cex=1.2)}
	mtext("Effect of Any Cue\n(95% Conf. Intervals)",side=1,line=3.5,cex=1.2)
	if(ommit.others==F){
			if(leg2){
			legend(x="bottom",
					legend=c("PSDBista","Other","Non-Partisan","PTista"),
					cex=1.2,pch=c(21,21,23,21),
					pt.bg=c(1,gray(0.5), gray(0.5),"white"),
					col=c(1,gray(0.5),gray(0.5),1),pt.cex=2,
					xpd=NA,inset=c(0,-0.23),bty="n",horiz=T)}
	}else{
			if(leg2){legend(x="bottom",
					legend=c("PSDBista","Non-Partisan","PTista"),
					cex=1.2,pch=c(21,23,21),
       				pt.bg=c(1,gray(0.5), "white"),
       				col=c(1,gray(0.5),1),pt.cex=2,
       				xpd=NA,inset=c(0,-0.23),bty="n",horiz=T) }
	}}
	if(port==T){
	if(leg){axis(side=2,at=1:5,tick=F,
			labels=c("5-Venezuela\nno Mercosul","4-Pré-sal",
					"3-Finan. Público\nEmp. Privadas","2-Proteção a\nInvestimento",
					"1-Salário\nMínimo"),las=2,cex=1.2)
		}
	mtext("Efeito de Algum Tratamento\n(Intervalos de Confiança de 95%)",side=1,line=3.5,cex=1.2)
	if(ommit.others==F){
		if(leg2){
		legend(x="bottom",
		legend=c("PSDBista","Outro","Sem Partido","PTista"),
		cex=1.2,pch=c(21,21,23,21),pt.bg=c(1,gray(0.5), gray(0.5),"white"),
		col=c(1,gray(0.5),gray(0.5),1),pt.cex=2,xpd=NA,
		inset=c(0,-0.23),bty="n",horiz=T)}
	}else{
	if(leg2){
	legend(x="bottom",
	legend=c("PSDBista","Sem Partido","PTista"),
	cex=1.2,pch=c(21,23,21),pt.bg=c(1,gray(0.5), "white"),
	col=c(1,gray(0.5),1),pt.cex=2,xpd=NA,inset=c(0,-0.23),bty="n",horiz=T) }
}}
}

plot.effects <- function(xx,leg=T,ommit.others=T,ommit.non=T,port=F){
	to.keep <- grep("45|13",colnames(xx$pes))
	pes <- xx[["pes"]][,to.keep]
	ses <- xx[["ses"]][,to.keep]
	offsets <- c(.25,.15,.05,-.05,-.15,-.25)
	pchs <- c(21,22,23)
	cols <- 1
	ltys <-1
	bgs <-  c(1,1,1,"white","white","white")	
	#if(sum(ommit.others,ommit.non)==2){
	#offsets <- c(.25,.15,.05,-.05,-.15,-.25)		
	#	}
	#if(sum(ommit.others,ommit.non)==1){
	#		
	#	}
	par(mar=c(9,8,.7,1))
	ylims <- c(0.5,0.5+nrow(pes))
	xlims <- c(-0.5,0.5)
	plot(xlims,ylims,type="n",ylab="",xlab="",yaxt="n")
	my.polygon<-function(x,col){#background
		polygon(x=c(xlims,xlims[2:1]),
				y=c(rep(x-0.5,2),rep(x+0.5,2)),
				col=gray(col),border=NA)
	}
	mapply(my.polygon,1:nrow(pes),
		rep(c(0.8, 0.7),nrow(pes))[1:nrow(pes)])
	for(i in 1:nrow(pes)){#loop over questions
		cilow <- pes[i,]-qnorm(.975)*ses[i,]
		cihigh <- pes[i,]+qnorm(.975)*ses[i,]
		segments(x0=cilow,x1=cihigh,y0=i+offsets,y1=i+offsets, #lines
				lwd=2,col=cols,lty=ltys)
		points(pes[i,],i+offsets,pch=pchs,bg=bgs,col=cols,cex=2) #points
		}		#end loop over questions
	abline(v=0,lty=2)
	if(port==F){
	axis(side=2,at=1:5,tick=F,
		labels=c("5-Venezuela\nin Mercosur","4-Offshore Oil\n(Pre-salt)",
				"3-Gvt. Financing\nPvt. Companies","2-Investment\nProtection",
				"1-Minimum\nWage"),las=2,cex=1.2)
	mtext("Effect of Each Treatment\n(95% Conf. Intervals)",side=1,line=3.3,cex=1.2)
	legend(x="bottom",legend=c("PSDBista X PTCue","PSDBista X PSDBCue",
				"PSDBista X DBLCue","PTista X PTcue","PTista X PSDBCue",
				"PTista X DBL Cue"),cex=1.2,pch=c(21,22,23,21,22,23),
				pt.bg=c(1,1,1,"white","white","white"),pt.cex=2,
				xpd=NA,inset=c(0,-0.35),bty="n",ncol=2)
	}
	if(port==T){
	axis(side=2,at=1:5,tick=F,
		labels=c("5-Venezuela\nno Mercosul","4-Pré-sal",
				"3-Finan. Público\nEmp. Privadas","2-Proteção a\nInvestimento",
				"1-Salário\nMínimo"),las=2,cex=1.2)
		mtext("Efeito de Cada Tratamento\n(Intervalos de Confiança de 95%)",side=1,line=3.3,cex=1.2)
		legend(x="bottom",legend=c("PSDBista X trat.\ PT",
			"PSDBista X trat.\ PSDB","PSDBista X ambos trat.","PTista X trat. PT",
			"PTista X trat PSDB","PTista X ambos trat."),cex=1.2,
			pch=c(21,22,23,21,22,23),pt.bg=c(1,1,1,"white","white","white"),
			pt.cex=2,xpd=NA,inset=c(0,-0.35),bty="n",ncol=2)
	}
}


plot.cue <- function(x,pty,cue="dbl",x2=NULL){
	#cue can be dbl,pt, or psdb
	#Takes output from estimate.pooled or estimate.cue
	#Plots barplots for partisans defines by pty, and treatement indicated by cue
	#Shows agreement with party whose object x was estimated
	#x2 is to add in both double and some cues in the same figure
	the.cue <- ifelse(cue=="dbl",":T4dblcue",
			   ifelse(cue=="pt",":T2ptcue",
			   ifelse(cue=="psdb",":T3psdbcue",
			   ifelse(cue=="some",":TT_somecue",break("Error in cue")))))
	baseline <- x$baseline[,paste("PP",pty,sep="")]
	the.order <- sort(baseline,index.return=T)$ix
	baseline <- baseline[the.order]
	effect <- x$pes[the.order,paste("PP",pty,the.cue,sep="")]
	tab <- rbind(baseline,baseline+effect)
	if(is.null(x2)==F){#if a x2 object is defined
		effect2 <- x2$pes[the.order,paste("PP",pty,":TT_somecue",sep="")]
		tab <- rbind(tab,baseline+effect2)
		}
	par(mar=c(4,4,1,1))
	xs <- barplot(tab,beside=T,names.arg = rep(NA,ncol(tab)),ylim=c(0,1)) 
	abline(h=seq(0.2,1,0.2),col=gray(0.8))
	barplot(tab,beside=T,names.arg = rep(NA,ncol(tab)),ylim=c(0,1),add=T) 
	ci.high <- baseline+effect+
				x$ses[the.order,paste("PP",pty,the.cue,sep="")]*qnorm(0.975)
	ci.low <- baseline+effect-
				x$ses[the.order,paste("PP",pty,the.cue,sep="")]*qnorm(0.975)
	segments(x0=xs[2,],x1=xs[2,], y0=ci.low, y1=ci.high,lwd=2)
	if(is.null(x2)==F){#if a x2 object is defined
		ci.high <- baseline+effect2+
				x2$ses[the.order,paste("PP",pty,":TT_somecue",sep="")]*qnorm(0.975)
		ci.low <- baseline+effect2-
				x2$ses[the.order,paste("PP",pty,":TT_somecue",sep="")]*qnorm(0.975)
		segments(x0=xs[3,],x1=xs[3,], y0=ci.low, y1=ci.high,lwd=2)
	}
	if(length(baseline)==5){
	my.labels=c("Venezuela\nin Mercosul","Offshore Oil",
			"Gvt. Financing\nPvt. Companies","Investment\nProtection",
			"Minimum\nWage")[the.order]
			}
	if(length(baseline)==6){
	my.labels=c("Venezuela\nin Mercosul","Offshore Oil",
			"Gvt. Financing\n(Online)","Investment\nProtection",
			"Minimum\nWage","Gvt. Financing...\n(BEPS)")[the.order]
			}	
	if(length(baseline)==1){my.labels=c("Gvt. Financing\nPvt. Companies\n(BEPS)")}
	axis(at=apply(xs,2,mean),side=1,labels=my.labels,line=1,tick=F,cex=0.75)
	ytext <- if(pty==13){"Share of PT identifiers agreeing with own party"}else{
				if(pty==45){"Share of PSDB identifiers agreeing with own party"}else{
				if(pty==0){"Share of non-partisans agreeing with PT"}else{
				if(pty==66){"Share of other partisans agreeing with with party"}	
				}}}
	mtext(ytext,side=2.5,line=3)
	mtext(paste('(with 95% conf. interval)'),side=2,line=2,cex=0.9)
	legend.cue <-ifelse(cue=="dbl","Double cue",
			   ifelse(cue=="pt","PT cue only",
			   ifelse(cue=="psdb","PSDB cue only",
			   ifelse(cue=="some","Some cue",break("Error in cue")))))
	if(is.null(x2)==T){cols <-gray.colors(2);
					   x2.cue <- NULL
					   }else{
					   cols <-gray.colors(3)
					   x2.cue <- "Any cue"}
	legend(x='topleft',legend=c("No cue",legend.cue,x2.cue),bty="n",	
			fill=cols,pt.cex=1.5)
}

plot.cue.H <- function(x,pty,cue="dbl",x2=NULL,cue2="some",sig=0.95){ #5% on one tailed test
	#Exactly the same as above, but plots results horizontally
	#Made small changes allowing for "second" bars to be anything, so this is more general now
	sig <- qnorm(sig)
	the.cue <- ifelse(cue=="dbl",":T4dblcue",
			   ifelse(cue=="pt",":T2ptcue",
			   ifelse(cue=="psdb",":T3psdbcue",
			   ifelse(cue=="some",":TT_somecue",break("Error in cue")))))
	the.cue2 <- ifelse(cue2=="dbl",":T4dblcue",
			   ifelse(cue2=="pt",":T2ptcue",
			   ifelse(cue2=="psdb",":T3psdbcue",
			   ifelse(cue2=="some",":TT_somecue",break("Error in cue")))))
	baseline <- x$baseline[,paste("PP",pty,sep="")]
	the.order <- sort(baseline,index.return=T)$ix
	baseline <- baseline[the.order]
	effect <- x$pes[the.order,paste("PP",pty,the.cue,sep="")]
	tab <- rbind(effect=baseline+effect,baseline)
	if(is.null(x2)==F){#if a x2 object is defined
		effect2 <- x2$pes[the.order,paste("PP",pty,the.cue2,sep="")]
		tab <- rbind(effect2=baseline+effect2,tab)
		}
	par(mar=c(5,7,0.5,1))
	xs <- barplot(tab,beside=T,names.arg = rep(NA,ncol(tab)),xlim=c(0,1.05),horiz=T) 
	rownames(xs) <- rownames(tab)
	abline(v=seq(0.2,1,0.2),col=gray(0.8))
	barplot(tab,beside=T,names.arg = rep(NA,ncol(tab)),xlim=c(0,1.05),horiz=T,add=T) 
	ci.high <- baseline+effect+
				x$ses[the.order,paste("PP",pty,the.cue,sep="")]*sig
	ci.low <- baseline+effect-
				x$ses[the.order,paste("PP",pty,the.cue,sep="")]*sig
	segments(y0=xs['effect',],y1=xs['effect',], x0=ci.low, x1=ci.high,lwd=2)
	if(is.null(x2)==F){#if a x2 object is defined
		ci.high <- baseline+effect2+
				x2$ses[the.order,paste("PP",pty,the.cue2,sep="")]*sig
		ci.low <- baseline+effect2-
				x2$ses[the.order,paste("PP",pty,the.cue2,sep="")]*sig
		segments(y0=xs['effect2',],y1=xs['effect2',], x0=ci.low, x1=ci.high,lwd=2)
	}
	if(length(baseline)==5){
	my.labels=c("Venezuela\nin Mercosul","Offshore Oil",
			"Gvt. Financing\nPvt. Companies","Investment\nProtection",
			"Minimum\nWage")[the.order]
			}
	if(length(baseline)==6){
	my.labels=c("Venezuela\nin Mercosul","Offshore Oil",
			"Gvt. Financing\nPvt. Companies\n(Online)","Investment\nProtection",
			"Minimum\nWage","Gvt. Financing\nPvt. Companies\n(BEPS)")[the.order]
			}	
	if(length(baseline)==1){my.labels=c("Gvt. Financing\nPvt. Companies\n(BEPS)")}
	axis(at=apply(xs,2,mean),side=2,labels=my.labels,line=0,tick=F,cex=1,las=2)
	ytext <- if(pty==13){"Share of PT identifiers agreeing with own party"}else{
				if(pty==45){"Share of PSDB identifiers agreeing with own party"}else{
				if(pty==0){"Share of non-partisans agreeing with selected party"}else{
				if(pty==66){"Share of other partisans agreeing with with party"}	
				}}}
	mtext(ytext,side=1,line=2.5,cex=1.2)
#	mtext(paste('(with 95% conf. interval)'),side=1,line=3.5,cex=0.9)
	legend.cue <-ifelse(cue=="dbl","Double cue",
			   ifelse(cue=="pt","PT cue only",
			   ifelse(cue=="psdb","PSDB cue only",
			   ifelse(cue=="some","Some cue",break("Error in cue")))))
	legend.cue2 <-ifelse(cue2=="dbl","Double cue",
			   ifelse(cue2=="pt","PT cue only",
			   ifelse(cue2=="psdb","PSDB cue only",
			   ifelse(cue2=="some","Some cue",break("Error in cue")))))
	if(is.null(x2)==T){cols <-gray.colors(2)[2:1];
					   x2.cue <- NULL
					   legend.cue2 <- NULL
					   }else{
					   cols <-gray.colors(3)[3:1]
					   x2.cue <- "Any cue"}
	legend(x='bottomright',legend=c("No cue",legend.cue,legend.cue2),bty="n",	
			fill=cols,pt.cex=1.5,cex=1.5)
}